% Estimate effects of US monetary policy shocks using model from
% Antolin-Diaz and Rubio-Ramirez (2018) under different combinations of
% (narrative) sign restrictions. Goal is to: 1) assess sensitivity of
% results to choice of conditional prior; and 2) disentangle
% informativeness of different restrictions.

clear variables
close all
clc

oldFolder = pwd;
cd ..
addpath([oldFolder,'/auxFunctions']);
cd(oldFolder);

%% Import data.
% Data are obtained from the replication files to Antolin-Diaz and
% Rubio-Ramirez (2018).
% The variables are:
% 1) Real GDP 2) GDP deflator 3) Commodity price index 4) Total reserves
% 5) Nonborrowed reserves 6) Federal funds rate
load('Uhlig_Data_Updated.mat');

% Re-order variables so that the shock of interest is ordered first. 
data = 100*data(:,[end 1:end-1]); % Federal funds rate ordered first

% Declare any exogenous variables (other than the constant; leave empty 
% if none).
exog = [];

% Set variable names for plotting figures.
varnames = {'Federal Funds Rate','Real GDP','GDP Deflator',...
    'Commodity Price Index','Total Reserves','Nonborrowed Reserves'};
% Set units for impulse responses. 0 = %, 1 = ppt, 2 = bps.
irfUnits = [2 0 0 0 0 0];

%% Options.
opt.p = 12; % No. of lags in VAR
opt.const = 1; % const = 1 if constant in VAR, = 0 otherwise
opt.jshock = 1; % Index of shock(s) of interest
opt.ivar = 2;  % Indices of variables of interest
opt.cumIR = []; % Indices of variables for which cumulative impulse response should be used
opt.H = 60; % Maximum horizon for impulse responses
opt.draws = 1000; % No. of draws from posterior of reduced-form parameters (phi) with non-empty identified set
opt.aalpha = 0.68; % Credibility level for credible intervals (0<aalpha<1)
opt.qMax = 100000; % Max no. of draws of orthonormal matrix Q when checking for empty identified set
opt.qDraws = 1000; % No. of draws of Q used to approximate bounds of (conditional) identified set
opt.dispIter = 100; % Print number of draws of phi remaining every dispIter draws
opt.gridLength = 1000; % Size of grid used when computing credible intervals

rng(23032021); % Set seed for random number generator

restr.elastRestr = []; % No elasticity restrictions

%% Sign restrictions only
% Each row of signRestr contains a vector (i,j,h,s,t) representing a
% `traditional' sign restriction, where t is the type of restriction:
% t = 1: the impulse response of the ith variable to the jth shock at the 
% hth horizon is nonnegative (s = 1) or nonpositive (s = -1).
% t = 2: the (ij)th element of A0 is nonnegative (s = 1) or nonpositive 
% (s = -1). 
% signRestr = []; % No sign restrictions on impulse responses
restr.signRestr = ...
      [1 1 0 1 1; % Response of FFR to monetary policy shock on impact is nonnegative
       1 1 1 1 1; % As above after one month
       1 1 2 1 1; % As above after two months
       1 1 3 1 1; % As above after three months
       1 1 4 1 1; % As above after four months
       1 1 5 1 1; % As above after five months
       3 1 0 -1 1; % Response of GDPDEF to monetary policy shock is nonpositive
       3 1 1 -1 1; % As above after one month, etc
       3 1 2 -1 1;
       3 1 3 -1 1;
       3 1 4 -1 1;
       3 1 5 -1 1;
       4 1 0 -1 1; % Response of COMM to monetary policy shock is nonpositive
       4 1 1 -1 1; % As above after one month, etc
       4 1 2 -1 1;
       4 1 3 -1 1;
       4 1 4 -1 1;
       4 1 5 -1 1;       
       6 1 0 -1 1; % Response of NBR to monetary policy shock is nonpositive
       6 1 1 -1 1; % As above after one month, etc
       6 1 2 -1 1;
       6 1 3 -1 1;
       6 1 4 -1 1;
       6 1 5 -1 1]; 

%% Run again with full set of sign and narrative restrictions.
% Each row of shockSignRestr contains a vector (i,s,t) representing the
% shock-sign restriction that the ith shock in period t is >= 0 (s = 1) or
% <= 0 (s = -1). t is equal to the index of the observation.
restr.shockSignRestr = ...
    [1, 1, find(dates == datenum(1974,4,01));
     1 1 find(dates == datenum(1979,10,01));
     1, 1, find(dates == datenum(1988,12,01));
     1, 1, find(dates == datenum(1994,2,01));
     1, -1, find(dates == datenum(1990,12,01));
     1, -1, find(dates == datenum(1998,10,01));
     1, -1, find(dates == datenum(2001,4,01));
     1, -1, find(dates == datenum(2002,11,01))];

% Each row of hdSignRestr contains a vector (i,j,t,h,k,s) representing a
% restriction on the historical decomposition. 
% k represents the class of restriction (Type A or Type B) and s represents
% the type of restriction within the class.
% k = 1: Type A restriction
%  - s = 1:  The jth shock is the 'most important contributor' to the 
%   change in the ith variable between periods t and  t+h.
%  - s = -1: The jth shock is the 'least important contributor' to the
%   change in the ith variable between periods t and t+h.
% k = 2: Type B restriction
%  - s = 1: The jth shock is the 'overwhelming contributor' to the change
%   in the ith variable between periods t and t+h.
% - s = -1: The jth shock is a 'negligible contributor' to the change in
%   the ith variable between periods t and t+h.
restr.hdSignRestr = ...
    [1, 1, find(dates == datenum(1974,4,01)), 0, 1, 1;
     1, 1, find(dates == datenum(1979,10,01)), 0, 1, 1;
     1, 1, find(dates == datenum(1988,12,01)), 0, 1, 1;
     1, 1, find(dates == datenum(1994,2,01)), 0, 1, 1;
     1, 1, find(dates == datenum(1990,12,01)), 0, 1, 1;
     1, 1, find(dates == datenum(1998,10,01)), 0, 1, 1;
     1, 1, find(dates == datenum(2001,4,01)), 0, 1, 1;
     1, 1, find(dates == datenum(2002,11,01)), 0, 1, 1];
% First shock in these periods was most important contributor to change in
% first variable (federal funds rate).

% Conduct (robust) Bayesian inference.
fprintf('Estimating model with sign restrictions, NR 6 and NR 7');

%% Construct data for estimating VAR for y_t.
YY = data(opt.p+1:end,:); % y_t 
XX = lags(data,1:opt.p); % Matrix of regressors in VAR for y_t
XX = XX(opt.p+1:end,:); % Drop initial missing observations
if opt.const == 1 % Add constant to matrix of regressors
    XX = [XX ones(size(XX,1),1)];
end
% Add exogenous variables to matrix of regressors
XX = [XX exog(opt.p+1:end,:)]; 

n = size(YY,2); % Number of variables
ni = length(opt.ivar); % Number of variables of interest
nj = length(opt.jshock); % Number of shocks of interest
m = size(XX,2); % Number of parameters in each equation for y_t
nExog = opt.const + size(exog,2); % Number of exogenous variables
T = length(YY); % Number of observations used in estimating VAR

% Adjust indices representing narrative restrictions to account for
% losing the first p observations. If the restrictions are imposed in the 
% first p periods of the original dataset, drop the restrictions.
if ~isempty(restr.shockSignRestr)
    restr.shockSignRestr(:,3) = restr.shockSignRestr(:,3) - opt.p;
    restr.shockSignRestr(restr.shockSignRestr(:,3) <= 0) = [];
end
if ~isempty(restr.hdSignRestr)
    restr.hdSignRestr(:,3) = restr.hdSignRestr(:,3) - opt.p;
    restr.hdSignRestr(restr.hdSignRestr(:,3) <= 0) = [];
end

%% Conduct posterior inference on impulse responses.
% j-th column of B gives MLE of coefficients in j-th reduced form VAR 
% equation for y_t. Last rows correspond to exogenous terms (if included).
phiHat.B = (XX'*XX)\XX'*YY;
phiHat.S = (YY - XX*phiHat.B)'*(YY-XX*phiHat.B);
phiHat.Sigma = (1/T)*phiHat.S; % MLE of innovation covariance matrix
phiHat.P = (XX'*XX)\eye(m);
phiHat.cholP = chol(phiHat.P,'lower');

% Storage arrays.
rMinPost = zeros(opt.draws,opt.H+1,ni,nj);
rMaxPost = zeros(opt.draws,opt.H+1,ni,nj);
rSinglePriorPost = zeros(opt.draws,opt.H+1,ni,nj);
B = zeros(n*m,opt.draws);
Sigma = zeros(n,n,opt.draws);

phiDraw = 0; % Counter for no. of draws from posterior of phi
nonStableCount = 0; % Counter for no. of draws with nonstable VAR
draw = 0; % Counter for no. of draws with non-empty IS
Qempty = [];
lbprob = zeros(opt.H+1);
ubprob = zeros(opt.H+1);

while draw < opt.draws
    
    % Posterior sampler with independent improper (Jeffreys) prior.
    phi.Sigma = iwishrnd(phiHat.S,T-m);
    phi.Sigmatr = chol(phi.Sigma,'lower');
    phi.Sigmatrinv = phi.Sigmatr\eye(n);
    phi.B = phiHat.B(:) + kron(phi.Sigmatr,phiHat.cholP)*randn(m*n,1);

    % Generate coefficients in orthogonal reduced-form VMA representation 
    % and check stability of VAR.
    [phi.vma,nonStable] = genVMA(phi,opt,nExog);
    nonStableCount = nonStableCount + nonStable;
    
    % If VAR is stable, proceed. Otherwise, redraw phi.
    
    if nonStable == 0
        
    phiDraw = phiDraw + 1;
    
    % Compute reduced-form VAR innovations.
    restr.U = (YY - XX*reshape(phi.B,m,n))';
        
    % Use Algorithm 1 (Step 2) to determine whether Q(phi|N,S) is nonempty
    % and draw Q satisfying restrictions.
    [Q0,restr.sphi,Qempty(phiDraw)] = drawQ0_ar(restr,phi,opt);
   
    if Qempty(phiDraw) == 0 % If Q(phi|N,S) is non-empty
            
            % Store reduced-form parameters.
            draw = draw + 1
            Sigma(:,:,draw) = phi.Sigma;
            B(:,draw) = phi.B;

            % Compute draw from posterior of impulse responses given 
            % conditionally uniform prior for Q|phi.
           
           % Find upper and lower bounds of conditional identified set for
           % impulse responses.
           [rMinPost,rMaxPost,etaDraw] = ...
               approximateBounds(restr,phi,opt);
           for hh=1:opt.H+1
               if rMinPost(hh)<0.0
                   ubprob(hh)=ubprob(hh)+1.0;
               end
               if rMaxPost(hh)<0.0
                   lbprob(hh)=lbprob(hh)+1.0;
               end 
           end
           % etaDraw(horizon,variable,shock,posterior draw)
           % Probabilities of decline for all Q: for output (2) and deflator (3) at horizons of 1, 2, 3, and 4 years

          
           if mod(opt.draws-draw,opt.dispIter) == 0
              
               fprintf('\n%d draws with non-empty identified set remaining...',...
               opt.draws-draw);
               
           end

    end
    
    end
    
end
ubprob = ubprob/opt.draws
lbprob = lbprob/opt.draws

% Table 7b:
disp('Table 7b')
disp([lbprob(13:12:61) ubprob(13:12:61)])

